Загрузим данные и пакеты

# Если у вас не установлены какие-то из библиотек ниже, то установить их можно следующей командой. Эти библиотеки пригодятся нам во время занятия, но подгружать в library() мы их будем тогда, когда будем затрагивать соответствующую тему.
# install.packages(c('ggpubr', 'plotly', 'corrplot', 'corrr', 'ggfortify', 'pheatmap', 'factoextra', 'FactoMineR', 'ggbiplot'))

# Загрузим библиотеки
library(dplyr)
## 
## Присоединяю пакет: 'dplyr'
## Следующие объекты скрыты от 'package:stats':
## 
##     filter, lag
## Следующие объекты скрыты от 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(ggpubr)

pima <- read.csv('pima.csv')

# Сделаем более детализированную переменную возрастных групп
pima <- pima %>% 
  mutate(
    diabetes_ch = as.character(diabetes),
    age_group = case_when(
      age < 31 ~ "21-30",
      age >= 31 & age < 41 ~ "31-40",
      age >= 41 & age < 51 ~ "41-50",
      age >= 51 & age < 61 ~ "51-60",
      age >= 61 ~ "60+"
    ))

table(pima$age_group)
## 
## 21-30 31-40 41-50 51-60   60+ 
##   417   157   113    54    27

…

Plotly - Интерактивные графики

Plotly - это специальная библиотека для создания интерактивных графиков со своим синтаксисом (очень похожим на ggplot). С ней можно познакомиться здесь.

Основные правила синтаксиса

Загрузим библиотеку:

# install.packages("plotly")
# или 
# devtools::install_github("ropensci/plotly")
library(plotly)
## 
## Присоединяю пакет: 'plotly'
## Следующий объект скрыт от 'package:ggplot2':
## 
##     last_plot
## Следующий объект скрыт от 'package:stats':
## 
##     filter
## Следующий объект скрыт от 'package:graphics':
## 
##     layout
skimr::skim(pima)
Data summary
Name pima
Number of rows 768
Number of columns 11
_______________________
Column type frequency:
character 3
numeric 8
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
diabetes 0 1 3 3 0 2 0
diabetes_ch 0 1 3 3 0 2 0
age_group 0 1 3 5 0 5 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
pregnant 0 1 3.85 3.37 0.00 1.00 3.00 6.00 17.00 ▇▃▂▁▁
glucose 0 1 120.89 31.97 0.00 99.00 117.00 140.25 199.00 ▁▁▇▆▂
pressure 0 1 69.11 19.36 0.00 62.00 72.00 80.00 122.00 ▁▁▇▇▁
triceps 0 1 20.54 15.95 0.00 0.00 23.00 32.00 99.00 ▇▇▂▁▁
insulin 0 1 79.80 115.24 0.00 0.00 30.50 127.25 846.00 ▇▁▁▁▁
mass 0 1 31.99 7.88 0.00 27.30 32.00 36.60 67.10 ▁▃▇▂▁
pedigree 0 1 0.47 0.33 0.08 0.24 0.37 0.63 2.42 ▇▃▁▁▁
age 0 1 33.24 11.76 21.00 24.00 29.00 41.00 81.00 ▇▃▁▁▁

Разберём базовые примеры. Для более подробного изучения приглашаю сюда

Box plot:

plot_ly(
  pima[pima$insulin != 0,],
  x = ~ insulin,
  color = ~ age_group,
  type = "box"
)

Scatter plot

plot_ly(data = pima[(pima$insulin != 0) & (pima$glucose != 0),],
        x = ~ insulin,
        y = ~ glucose)
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Начнём усложнять

Добавим название и настройки маркеров:

plot_ly(
  data = pima[(pima$insulin != 0) & (pima$glucose != 0),],
  x = ~ insulin,
  y = ~ glucose,
  marker = list(
    size = 10,
    # размер
    color = 'rgba(255, 182, 193, .9)',
    # Цвет внутри
    line = list(color = 'rgba(152, 0, 0, .8)',  # Цвет окружности
                width = 2)
  )
)   %>%
  layout( #настраиваем вывод подписей
    title = 'Отношение уровня глюкозы и инсулина в данных PIMA',
    yaxis = list(title = 'Уровень глюкозы',
                 zeroline = FALSE),  # Уберём выделения нулевых осей по y
    xaxis = list(title = 'Уровень инсулина',
                 zeroline = FALSE)) # Уберём выделения нулевых осей по y
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Раскрасим по диабет-статусу:

plot_ly(
  data = pima[(pima$insulin != 0) & (pima$glucose != 0),],
  x = ~ insulin,
  y = ~ glucose,
  color = ~diabetes_ch
)   %>%
  layout(
    title = 'Отношение уровня глюкозы и инсулина в данных PIMA',
    yaxis = list(title = 'Уровень глюкозы',
                 zeroline = FALSE),  # Уберём выделения нулевых осей по y
    xaxis = list(title = 'Уровень инсулина',
                 zeroline = FALSE)) # Уберём выделения нулевых осей по y
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

Добавим третью ось:

plot_ly(
  data = pima[(pima$insulin != 0) & (pima$glucose != 0) & (pima$mass != 0),],
  x = ~ insulin,
  y = ~ glucose,
  z = ~mass,
  color = ~diabetes_ch
)   %>%
  layout(
    title = 'Отношение уровня глюкозы и инсулина в данных PIMA',
    yaxis = list(title = 'Уровень глюкозы',
                 zeroline = FALSE),  # Уберём выделения нулевых осей по y
    xaxis = list(title = 'Уровень инсулина',
                 zeroline = FALSE)) # Уберём выделения нулевых осей по y
## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

Palette / Цветовые гаммы

# Со стандартной палеткой
plot_ly(data = pima[(pima$insulin != 0) & (pima$glucose != 0),],
        x = ~ insulin,
        y = ~ glucose,
        color = ~age_group)
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
# С изменённой палеткой
plot_ly(data = pima[(pima$insulin != 0) & (pima$glucose != 0),],
        x = ~ insulin,
        y = ~ glucose,
        color = ~age_group,
        colors = "Set1")
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Палетки можно сделать самостоятельно (аналогичным образом это делается в ggplot):

pal <- c("red", "blue", "green")
pal <- setNames(pal, levels(pima$age_group)) #age_group - должнабыть факторная. Число цветов в pal должно соответствовать числу уровней фактора.

plot_ly(data = pima[(pima$insulin != 0) & (pima$glucose != 0),],
        x = ~ insulin,
        y = ~ glucose,
        color = ~age_group,
        colors = pal)
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Цвет, размер точек и подписи к ним

plot_ly(
  data = pima[(pima$insulin != 0) & (pima$glucose != 0),], 
  x = ~insulin, 
  y = ~glucose,
  color = ~diabetes_ch, 
  size = ~mass,
  text = ~mass, 
  hoverinfo = "text")
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

Анимированный график

Построим анимацию для разных возрастов. Она не будет особо осмысленна (т.к. делается обычно для данных по времени), но мы учимся:

plot_ly(
  data = pima[(pima$insulin != 0) & (pima$glucose != 0),],
  x = ~insulin, 
  y = ~glucose, 
  color = ~diabetes_ch, 
  size = ~mass,
  text = ~mass, 
  frame = ~age, # По какой колонке нам сделать анимацию/интерактивность?
  hoverinfo = "text",
  type = 'scatter',
  mode = 'markers'
) 
## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.

## Warning: `line.width` does not currently support multiple values.
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

Почему мы не учим plotly в первую очередь?

Он не гибкий. Если ggplot понимает что вам нужно на ходу, то для plotly вам нужно подготовить данные в том виде в котором вы хотите их визуализировать.

ЗАДАНИЕ

В чанке ниже напишите код 3d графика ао осям age, mass, pressure:

plot_ly(
  data = pima[(pima$age != 0) & (pima$mass != 0) & (pima$pressure != 0),],
  x = ~ age,
  y = ~ mass,
  color = ~pressure
)
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Обход - ggplotly

Когда вы хотите представить итоговый результат в интерактивном формате (например, с помощью странички на сайте сделанной в Rmd), то её лучше писать в plotly. Однако, если вы работаете в формате эксплораторного анализа, то можно воспользоваться “обёрткой” plotly для вашего ggplot графика. Сделать это очень просто.

Сделаем и сохраним обычный график:

plot <- pima %>% 
  filter(mass != 0 & triceps != 0) %>% 
  ggplot(aes(x=mass, y=triceps, color = diabetes)) + 
  geom_point(size=3) +
  theme_minimal()

“Обернём” наш график:

ggplotly(plot)

ЗАДАНИЕ

В чанке ниже напишите код любого графика и сделайте его интерактивным:

insurance_cost <- read.csv('../hw1/data/insurance_cost.csv', stringsAsFactors = T)
plot1 <- insurance_cost %>% 
  #filter(mass != 0 & triceps != 0) %>% 
  ggplot(aes(
    x = age,
    y = charges,
    color = smoker,
    fill = smoker,
    group = smoker)) + 
    geom_point(size=3) +
  geom_smooth(method=lm, # Выбираем модель. Стандартно - линейная модель
              color="red", fullrange = T,
              fill="#69b3a2", 
              se=TRUE) + # Используем ли доверительные интервалы?
 theme_minimal() +
 ggtitle('Scatter plot отношения переменных возраста и зарплаты') +
 theme(axis.text.x = element_text(size = 14))

ggplotly(plot1)
## `geom_smooth()` using formula = 'y ~ x'

…

Анализ корреляций

Когда мы говорим об отношениях переменных друг к другу, мы так или иначе приходим к языку корреляции. Т.е. статистической взаимосвязи двух или более случайных величин при которых изменения в одной переменной связаны с изменением во второй.

На всякий случай привожу формулу рассчёта коэффициента корреляции Пирсона: \(r_{xy}=\frac{\Sigma(x_i-\bar{x})\times(y_i-\bar{y})}{\sqrt{\Sigma(x_i-\bar{x})^2\times\Sigma(y_i-\bar{y})^2}}\)

Он может принимать значения от -1 (отрицательная связь), до 1 (положительная связь).

В R можно строить матрицу корреляций для всех численных переменных нашего датасета с помощью пакета corrplot. Загрузим его:

library(corrplot)
## corrplot 0.92 loaded

Сама матрица строится в 2 этапа.

  1. Получаем объект матрицы:
# Для более "чистого" результата, избавляемся от ошибочных значений
pima_clear <- pima %>% 
  filter(glucose != 0 & pressure != 0 & triceps != 0 & insulin != 0 & mass != 0 & age != 0 ) %>% 
  select(is.integer | is.numeric) # Обратите внимание, в dplyr можно задавать выборку колонок через команды определения формата данных
## Warning: Use of bare predicate functions was deprecated in tidyselect 1.1.0.
## ℹ Please use wrap predicates in `where()` instead.
##   # Was:
##   data %>% select(is.integer)
## 
##   # Now:
##   data %>% select(where(is.integer))
## Warning: Use of bare predicate functions was deprecated in tidyselect 1.1.0.
## ℹ Please use wrap predicates in `where()` instead.
##   # Was:
##   data %>% select(is.numeric)
## 
##   # Now:
##   data %>% select(where(is.numeric))
head(pima_clear)
# Получаем непосредственно матрицу
pima_cor <- cor(pima_clear)
pima_cor
##              pregnant   glucose   pressure   triceps    insulin        age
## pregnant  1.000000000 0.1982910  0.2133548 0.0932094 0.07898363 0.67960847
## glucose   0.198291043 1.0000000  0.2100266 0.1988558 0.58122301 0.34364150
## pressure  0.213354775 0.2100266  1.0000000 0.2325712 0.09851150 0.30003895
## triceps   0.093209397 0.1988558  0.2325712 1.0000000 0.18219906 0.16776114
## insulin   0.078983625 0.5812230  0.0985115 0.1821991 1.00000000 0.21708199
## age       0.679608470 0.3436415  0.3000389 0.1677611 0.21708199 1.00000000
## mass     -0.025347276 0.2095159  0.3044034 0.6643549 0.22639652 0.06981380
## pedigree  0.007562116 0.1401802 -0.0159711 0.1604985 0.13590578 0.08502911
##                 mass     pedigree
## pregnant -0.02534728  0.007562116
## glucose   0.20951592  0.140180180
## pressure  0.30440337 -0.015971104
## triceps   0.66435487  0.160498526
## insulin   0.22639652  0.135905781
## age       0.06981380  0.085029106
## mass      1.00000000  0.158771043
## pedigree  0.15877104  1.000000000
  1. Визуализируем её в corplot:
corrplot(pima_cor, method = 'number')

У представления такой матрицы может быть ещё большое количество вариаций. Их можно посомтреть здесь.

Кроме того, я рекомендую посмотреть также на пакет corrr.

В нём реализованы более дателизированные функции для анализа корреляций. Например функция ниже показывает не только корреляционные взаимосвязи, но и отношения близости переменных друг к другу с точки зрения сетевого анализа:

library(corrr)

pima_cor %>% 
  network_plot(min_cor = .0)

…

Heat map

Анализ двух категориальных переменных

В базовом варианте, heat map используют для иллюстрации отношения двух категориальных переменных друг с другом. Для этого, нужно сначала агрегировать их значения:

pima_aggr <- pima %>% 
  group_by(age_group, diabetes) %>% 
  summarise(N = n())
## `summarise()` has grouped output by 'age_group'. You can override using the
## `.groups` argument.

Визуализируем:

pima_aggr %>% 
  ggplot(aes(x = age_group, y = diabetes, fill = N)) +
  geom_tile(color = "black") +
  geom_text(aes(label = N), color = "white", size = 4) +
  coord_fixed()

Анализ наблюдений

Второй способ использовать heat map - представлять сами наблюдения по какому-то универсальному признаку.

Например, мы можем применить процедуру стандартизации значений для всех наблюдений по численным переменным и визуализировать его.

Стандартизируем значения:

pima_clear_scaled <- scale(pima_clear)
head(pima_clear_scaled)
##         pregnant    glucose    pressure    triceps    insulin        age
## [1,] -0.71651083 -1.0896533 -0.37317791 -0.5843629 -0.5221747 -0.9670632
## [2,] -1.02789913  0.4657189 -2.45382847  0.5567094  0.1005024  0.2093178
## [3,] -0.09373423 -1.4460927 -1.65357826  0.2714413 -0.5726620 -0.4769045
## [4,] -0.40512253  2.4099341 -0.05307782  1.5076030  3.2559608  2.1699528
## [5,] -0.71651083  2.1507054 -0.85332804 -0.5843629  5.8055711  2.7581434
## [6,]  0.52904237  1.4054229  0.10697222 -0.9647204  0.1594043  1.9738893
##            mass   pedigree
## [1,] -0.7095143 -1.0305593
## [2,]  1.4249091  5.1085822
## [3,] -0.2968591 -0.7961084
## [4,] -0.3680065 -1.0566094
## [5,] -0.4249245 -0.3619399
## [6,] -1.0367925  0.1851123

Визуализируем с помощью расширения к ggplot - ggfortify(). В чем проблема такого графика?

library(ggfortify) 

autoplot(pima_clear_scaled)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

…

Методы уменьшения размерности

Поиск схожих наблюдений. Лирическое отступление об основах кластерного анализа

Иерархическая кластеризация (hierarchical clustering)

Относится к семейству аггломеративных кластеризаций (agglomerative/agnes, есть ещё другой вид - divisive/diana), но является самой простой из них - не требует выделить число кластеров заранее. Нужна в первую очередь, чтобы увидеть сруктуру отнесения наблюдений к кластерам (и подумать о том, как вообще устроенны наши наблюдения).

Принцип работы - каждое наблюдение находит ближайшее себе по дистанции до него (по каждой из переменных). После нахождения первого соседа, начинатся следующий шаг, когда образовавшийся мини-кластер ищет другие ближайшие к ним кластеры/наблюдения. Так прохожит несколько шагов, пока все наблюдения не сольются в один кластер.

“Длинна” дистанции обозначает меру схожести наблюдений - Высоту/Height (Площадь Восстания ближе к ЕУ чем Кудрово).

Посмотрим визуально на картинке.

…

Есть много способов измерять расстояние. Мы будем пользоваться базовым - Евклидовым расстоянием.

Как с этим работать?

  1. Установик и загрузим библиотеку для визуализации кластеров
library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
  1. Создаём матрицу дистанций
pima_clear_dist <- dist(pima_clear_scaled, method = "euclidean")
as.matrix(pima_clear_dist)[1:6,1:6]
##          1        2        3        4        5        6
## 1 0.000000 7.222446 1.831606 6.407444 8.073168 4.341744
## 2 7.222446 0.000000 6.633504 8.115465 8.899023 6.752575
## 3 1.831606 6.633504 0.000000 6.386919 8.127109 4.604262
## 4 6.407444 8.115465 6.386919 0.000000 3.537396 4.431248
## 5 8.073168 8.899023 8.127109 3.537396 0.000000 6.028399
## 6 4.341744 6.752575 4.604262 4.431248 6.028399 0.000000
  1. Высчитываем дендрограмму кластеров. Обратите внимание на аргумент method. Методов для кластеризации достаточно много. Мы используем один из базовых - “ward.D2”
pima_clear_hc <- hclust(d = pima_clear_dist, 
                        method = "ward.D2")
  1. Визуализируем (на больших данных может занять какое-то время)
fviz_dend(pima_clear_hc, 
          cex = 0.1) # cex() - размер лейблов
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <]8;;https://github.com/kassambara/factoextra/issueshttps://github.com/kassambara/factoextra/issues]8;;>.

Далее можно углубляться в способы визуализации и оценки таких результатов. Рекомендую посмотреть на учебник по ссылке.

…

Heat map + Tree map?

Мы затронули основы кластеризации, чтобы понять, как работать в тех случаях, когда нам нужно упорядочить данные, сгруппировав наблюдения. Совместив это с heat map мы можем выделить группы наблюдений по каким-либо признакам.

Сделать это просто через пакет pheatmap.

Загрузим его.

library(pheatmap)

Визуализируем heat map и tree map одновременно:

pheatmap(pima_clear_scaled)

Что не так в этом графике?

Principal component analysis (PCA)

Сейчас мы уже можем выделить определённые группы пациентов. Однако даже на таком маленьком объёме данных (напомню, мы работаем с 392 наблюдениями) график почти не читаем. Что делать?

Одно из самых популярных решений - метод главных компонент (principal component analysis или PCA).

…

Теория и учебный пример

Чтобы разбраться с методом рассмотрим какие две проблемы он решает и как он это делает:

  1. Если в данных есть ряд скоррелированных переменных (например, в нашем случае - mass и triceps), то при анализе они не дадут провести корректные оценки данных и выделение паттернов.
Почему?
  1. Если в данных очень много наблюдений и переменных, то при разведывательном анализе необходимо сократить их размер, не потеряв в качестве и полноте интерпретаций. Эту задачу выполняют методы уменьшения размерности и PCA относится к ним.

…

Рассмотрим учебный пример. Допустим у нас в данных есть только колонки mass и triceps. Мы хотим с одной стороны снять корреляцию между ними, а с другой - уменьшить размерность, т.е. представить их одной переменной.

Как нам это сделать?

…

Мы можем сконструировать новые переменные из двух имеющихся и сняв корреляцию и уменьшив размерность!

Давайте сначала ещё раз посмотрим на данные:

# Подготовим данные
pima_example <- pima_clear %>% 
  select(mass, triceps)

# Визуализируем
ggplot() +
  geom_point(data = pima_example, aes(x = mass, y = triceps)) +
  theme_minimal()

# Что значит такое распределение точек?

Сконструируем переменные, в которых не будет корреляции и которые (делаем из бага фичу!) помогут уменьшить размерность

Как нам сделать из двух переменных одну?

Если мы хотим привести две переменные к одной, то мы должны провести какое-то действие с этими переменными. Например, сложить. Давайте сделаем это:

pima_example <- pima_example %>% 
  mutate(pc1 = mass + triceps) # Создаем первую главную компоненту

Давайте посмотрим, как на данных выглядят три типа значений. Те, которые близки к медиане по pc1, те которые близки к максимуму по pc1 и те, которые близки к минимуму по pc1.

Создадим дополнительную колонку по этим значениям:

pima_example <- pima_example %>% 
  mutate(pc1_decile = ntile(pc1, 10))

А теперь визуализируем наш scatter plot соотношения mass и triceps, НО, раскрасим цветом значения первого, пятого, шестого и десятого децелей:

Подготовим данные, добавив в них номинативную переменную исходных децелей:

pima_example <- pima_example %>% 
  mutate(pc1_decile_ch = case_when(
    pc1_decile == 1 ~ "1",
    (pc1_decile == 5) | (pc1_decile == 6) ~ "5-6",
    pc1_decile == 10 ~ "10"
  ))
pima_example
ggplot() +
  geom_point(data = pima_example, 
             aes(x = mass, 
                 y = triceps, 
                 color = pc1_decile_ch)) +
  theme_minimal()

# Что мы видим?

Такой результат мы получаем на сложении. А давайте теперь поэксперементируем с вычитанием:

pima_example <- pima_example %>% 
  mutate(pc2 = mass - triceps) %>% # Создаём колонку для результатов вычитания
  mutate(pc2_decile = ntile(pc2, 10)) %>% # Находим её децили
    mutate(pc2_decile_ch = case_when( # Отмечаем номинативной переменной значения разных децелей
    pc2_decile == 1 ~ "1",
    (pc2_decile == 5) | (pc2_decile == 6) ~ "5-6",
    pc2_decile == 10 ~ "10"
  ))

Визуализируем:

ggplot() +
  geom_point(data = pima_example, 
             aes(x = mass, 
                 y = triceps, 
                 color = pc2_decile_ch)) +
  theme_minimal()

# Что мы видим?

А теперь завершим нашу логику. Давайте оставим в данных только средние значения по каждому квантилю и потом нарисуем по ним линию тренда.

Создадим данные:

pima_example_pc1 <- pima_example %>%
  group_by(pc1_decile) %>% # Группируем по квантилям (1:10)
  summarise(mass_pc1 = mean(mass), # В каждом квантиле находим среднее для mass
            triceps_pc1 = mean(triceps)) # В каждом квантиле находим среднее для triceps

# Сделаем тоже самое с pc2
pima_example_pc2 <- pima_example %>%
  group_by(pc2_decile) %>%
  summarise(mass_pc2 = mean(mass),
            triceps_pc2 = mean(triceps))

Теперь визуализируем график по mass и triceps, НО добавив линию тренда по средним значениям pc1 и pc2 для mass и triceps:

ggplot() +
  geom_point(data = pima_example,
             aes(x = mass, y = triceps)) + # Основные данные
  geom_smooth(data = pima_example_pc1, # Отмечаем линию тренда по pc1
              aes(x = mass_pc1, y = triceps_pc1), # Используем значения mass и triceps, которые являются средними для децилей pc1
              method=lm,
              color="Blue", fullrange = F,
              size = 2
              ) +
  geom_smooth(data = pima_example_pc2, # Аналогично с pc2
              aes(x = mass_pc2, y = triceps_pc2),
              method=lm, 
              orientation = "y", # Технический трюк. Просим R отстраивать линию вдоль оси y, а не x
              color="Green", fullrange = F,
              size = 2
              ) +
  theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

А теперь, следим за руками...

А что если мы просто сделаем эти линии тренда нашими новыми координатными осями?
Т.е. вместо x будет pc1, вместо y - pc2, а точки перерисуются исходя из нашего нового представления.

Вспоминаем две проблемы, с которых мы начинали. Как такое преобразование помогает решить обе сложности?

…

Чтобы не быть голословными, перерисуем график, но уже с новыми координатами - pc1 и pc2:

ggplot() +
  geom_point(data = pima_example, aes(x = pc1, y = pc2)) +
  theme_minimal() 

# Догодайтесь, где теперь лежат средние линии наших старых координат (mass, triceps)?
Вопросы на понимание:

* А зачем мы это сделали? Какие следствия такого подхода?

* Если мы хотим выразить две наши переменные какой-то одной, то что мы выберем pc1 или pc2?

…

Подойдём к вопросу более формально

  • Нам нужно найти оси, которые приведут корреляцию к нулю и при этом, создадут макисмальную полноту разброса переменных. Т.е. фактически, в нашей задаче мы решем два уравнения:

$ PC1 = (alpha_1 * mass) + (beta_1 * triceps) $ $ PC2 = (alpha_2 * mass) + (beta_2 * triceps) $

где alpha и beta - это коэффициенты, которые позволяют добиться минимальной скоррелированности mass и triceps. (читай - найти такой угол для новых осей координат, где в данных будет максимальный разброс при минимальной корреляции)

Воспользуемся специальной функцией для этого:

pima_example <- pima_example %>% 
  select(mass, triceps) # Оставим только нужные нам две переменные

pima.pca <- prcomp(pima_example, 
                scale = T) # Нужно ли нормировать? TRUE|FALSE. Фактически, мы могли бы подать на функцию pima_example_scaled и тогда следовало поставить scale = F, веть нормирование в том датафрейме уже сделано!

Посмотрим на веса наших главных компонент:

pima.pca$rotation
##               PC1        PC2
## mass    0.7071068 -0.7071068
## triceps 0.7071068  0.7071068

Rotation показывает коэффициенты разворота новых осей относительно старых. Таким образом, решением нашего уравнения будет:

$ PC1 = (0.707 * mass) + (-0.707 * triceps) $ $ PC2 = (0.707 * mass) + (0.707 * triceps) $

(странность весов обусловлена тем, что мы работаем с учебным примером, где всего лишь 2 переменные)

Нарисуем компоненты на графике:

ggplot() +
  geom_point(data = as.data.frame(pima.pca$x), # Новые переменные лежат здесь
             aes(x = PC1, y = PC2)) +
  theme_minimal()

Но как функция высчитала эти цифры? 1. Она находит центр распределения данных по среднему двух переменных; 2. От него строит множество прямых линий с каждым разом делая угол наклона всё больше; 3. Для каждой линии считается сумма квадратов расстояний от всех точек в данных до линии; 4. Выбирается линия, где эта сумма будет минимальной. Она и становится новой главной компонентой; 5. Если алгоритм уже отобрал однй главную компоненту, он начинает отбирать следующую так, чтобы каждая следующая линия имела наименьшую корреляцию со всеми предыдущими; 6. Алгоритм останавливается, тогда, когда количесто найденных главных компонент не сравняется с количеством переменных.

Но в чём суть, если мы хотели уменьшить количество переменных, а оно (количество главных компонент) осталось таким же?

- Да, осталось, но только вот главные компоненты концептуально неравны друг другу. Каждая следующая из них объясняет меньше предыдущей. Считается, что алгоритм сработал хорошо, когда 70% вариации данных укладывается в 3 компоненты. 

...Т.е. компонент может быть десятки, но первые 3 агреггируют в себе большую часть сложности наших данных.

Это легко заметить уже даже в случае трёх переменных (и трёх компонент). Рассмотрим прекрасный материал Гарика Мороза и соавторов по ссылке здесь.


```r
# Кстати, мы теперь также можем делать 3d визуализацию с помощью plotly:
plot_ly(data = pima_clear, 
        x=~mass, 
        y=~pressure, 
        z=~triceps, 
        size = 1,
        type="scatter3d", mode="markers")
```

```{=html}
<div id="htmlwidget-e7fe1dd198e710c2c86a" style="width:672px;height:480px;" class="plotly html-widget"></div>
<script type="application/json" data-for="htmlwidget-e7fe1dd198e710c2c86a">{"x":{"visdat":{"2810371366e6":["function () ","plotlyVisDat"]},"cur_data":"2810371366e6","attrs":{"2810371366e6":{"x":{},"y":{},"z":{},"mode":"markers","size":1,"alpha_stroke":1,"sizes":[10,100],"spans":[1,20],"type":"scatter3d"}},"layout":{"margin":{"b":40,"l":60,"t":25,"r":10},"scene":{"xaxis":{"title":"mass"},"yaxis":{"title":"pressure"},"zaxis":{"title":"triceps"}},"hovermode":"closest","showlegend":false},"source":"A","config":{"modeBarButtonsToAdd":["hoverclosest","hovercompare"],"showSendToCloud":false},"data":[{"x":[28.1,43.1,31,30.5,30.1,25.8,45.8,43.3,34.6,39.3,36.6,31.1,23.2,22.2,31.6,24.8,24,37.1,34,45.4,19.4,24.2,24.4,33.7,34.7,37.7,46.8,41.5,25.4,19.6,28.9,32.9,28.6,35.1,29.3,32.4,38.5,37.1,32,46.7,24.7,33.9,20.4,28.7,49.7,26.6,28.7,29.5,34.3,37.4,33.3,34,31.2,30.5,23.2,53.2,33.6,55,42.9,33.3,34.5,29.7,34.5,21.1,33.8,30.8,28.7,36.9,32.5,32.8,30.5,37.4,34.3,40.6,24.6,25.2,29,40.9,37.2,44.2,29.9,28.4,35.4,43.5,29.7,32.7,67.1,34.9,30.1,32,27.9,31.6,33.1,39.4,22.9,34.8,30.9,20.4,37.7,37.5,33.2,42.6,34.2,41.8,35.8,30,34.6,28.8,23.6,34.6,36.7,45.2,46.2,25.4,29.7,35.9,33.1,27.1,38.2,52.3,35.4,24.4,27.6,25.9,33.3,30.9,33.6,25.2,33.2,40.5,27.8,25.3,35.9,32.4,26,38.7,45.6,20.8,36.1,36.9,36.6,43.3,40.5,35.5,28,30.7,36.6,31.6,35.8,39.7,25.5,24.8,30.5,32.9,39.4,26.6,29.5,34.1,19.3,38.1,27.5,26.8,25.6,35.1,45.5,30.8,32.7,23.9,47.9,34.2,25.9,25.9,38.5,28.7,21.8,27.2,33.3,35.3,36.5,31.2,34.9,34,27.5,32.8,38.4,35.8,34.9,36.2,39.2,25.2,37.2,43.4,30.8,25.4,25.1,24.3,22.3,32,31.6,32,23.7,22.1,27.7,24.7,35,42.1,42.4,34.4,42.4,26.2,34.6,35.7,26.4,45.3,26,40.6,42.9,37,34.1,40.6,35,30.4,30,32.2,33.2,59.4,25.3,36.5,33.6,30.5,21.2,39.9,37.8,30.2,37.6,25.9,20.8,35.3,21.8,27.8,36.8,46.1,33.7,23.8,25.9,35.5,27.8,38.2,42.3,40.7,46.5,36.8,28.9,30.1,25.1,29.3,25.2,33.3,36.5,28.6,30.4,22.1,25.6,31.6,30.3,19.6,25,33.2,18.2,26.3,30.8,29.8,41.3,33.3,36.3,36.4,39.4,32.4,39.5,32,34.5,43.6,33.1,32.8,31.9,29.9,36.9,25.5,41.3,37.6,26.9,26.1,38.6,32,31.3,34.3,29.5,34.7,30.1,35.5,24,28.7,33.3,39.4,28.5,33.6,32,27.8,23.1,35.2,40,19.5,41.5,24,30.9,32.9,38.2,36.1,20.1,38.4,43.5,37.7,34.5,27.5,31.6,40.9,19.5,29.3,27.6,39.4,23.4,37.8,28.3,25.2,33.8,34.1,34.2,38.7,21.8,38.9,34.2,37.6,37.9,34.8,34,30.9,33.6,35.5,57.3,24.2,24.2,44.6,33.2,24.1,46.1,39.1,38.5,30.4,29.9,34.5,35.9,28.4,34.4,38,31.2,29.6,26.4,33.9,33.8,35.5,38.1,29.3,39.1,36.1,28.4,44.5,29,27.4,36.6,42.3,30.8,28.5,40.6,30,46.3,36.4,39,43.3,36.5,28.4,32.9,26.2],"y":[66,40,50,70,60,72,84,30,70,88,94,70,66,82,76,58,60,72,64,110,80,50,66,90,66,68,88,64,58,66,85,66,64,86,78,74,68,70,80,78,82,72,48,50,90,72,56,58,58,85,72,62,76,54,76,76,74,30,70,58,88,70,64,68,60,70,60,72,52,62,64,74,86,82,52,56,74,72,74,80,74,90,70,60,64,72,110,64,68,98,76,80,70,84,62,64,60,70,72,76,64,65,82,70,62,68,60,60,66,78,70,80,80,80,68,84,70,50,76,90,70,80,62,50,76,68,74,62,78,70,64,62,76,88,74,84,86,56,72,88,62,78,48,62,70,84,78,58,82,76,68,68,68,68,70,74,50,68,80,66,60,90,72,64,86,70,58,60,76,78,70,74,88,46,62,62,50,74,76,64,74,54,86,102,82,64,58,52,82,82,60,100,72,60,62,70,54,82,68,66,64,72,58,56,84,48,68,72,84,74,60,84,64,88,68,64,78,78,64,94,82,74,74,66,64,78,72,80,64,74,64,68,54,68,84,74,72,70,56,52,64,78,80,76,74,70,58,82,68,62,78,65,70,72,70,74,90,64,90,60,50,62,54,70,88,90,70,80,64,74,66,60,66,56,80,92,74,72,90,78,90,76,68,82,68,62,64,70,66,68,60,54,72,62,72,66,58,60,86,44,44,76,86,78,52,72,82,24,38,78,78,62,82,62,54,58,88,74,62,86,70,88,78,82,76,76,74,86,72,74,74,50,84,54,60,74,70,52,58,80,82,106,80,80,58,78,68,106,100,58,56,64,74,74,82,70,68,90,74,88,76,76,46,64,64,78,58,50,78,60,66,68,86,78,78,88,56,86,60,80,44,58,88,84,74,70,78,88,88,58,76,72],"z":[23,35,32,45,23,19,47,38,30,41,33,26,15,19,36,11,33,47,25,24,11,15,21,34,42,39,60,41,34,13,27,20,35,20,26,29,25,32,15,40,18,27,18,30,51,18,29,28,31,25,33,26,34,32,15,56,30,42,30,36,24,14,37,13,20,26,25,29,26,31,35,50,28,42,15,21,19,41,40,34,18,12,23,42,24,42,46,18,36,41,39,35,44,41,13,44,27,16,32,29,27,26,32,40,41,30,29,33,15,27,39,31,37,25,28,21,32,22,35,33,33,14,7,16,28,15,18,32,50,52,23,10,28,15,26,44,39,17,43,29,30,37,45,31,38,29,25,33,41,37,23,14,19,28,37,17,10,22,11,39,12,33,21,32,36,32,16,18,43,34,13,21,36,19,19,12,40,40,36,33,25,28,16,28,48,22,40,43,43,15,37,39,30,8,18,24,13,26,23,29,14,12,24,34,41,32,49,30,23,22,35,33,29,41,18,46,32,39,30,46,25,16,11,23,27,63,12,45,37,18,13,32,28,28,48,33,22,40,13,10,36,41,45,17,38,30,22,31,42,41,32,28,18,15,33,32,19,25,26,23,23,17,19,18,34,7,32,33,19,15,31,18,52,30,37,49,40,25,23,29,35,27,21,43,30,24,23,33,32,34,19,14,30,32,29,30,31,17,30,47,20,24,27,50,22,45,14,19,18,29,42,25,39,13,21,22,42,26,13,42,27,47,40,17,18,32,12,17,30,35,17,36,35,25,23,40,28,27,35,48,31,46,46,45,33,30,26,23,35,17,28,39,26,26,46,32,49,24,19,11,27,20,21,32,13,27,20,33,39,46,36,29,30,29,23,37,27,27,17,37,20,18,37,33,41,22,39,44,39,26,48,23],"mode":"markers","type":"scatter3d","marker":{"color":"rgba(31,119,180,1)","size":[55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55,55],"sizemode":"area","line":{"color":"rgba(31,119,180,1)"}},"textfont":{"size":55},"error_y":{"color":"rgba(31,119,180,1)","width":55},"error_x":{"color":"rgba(31,119,180,1)","width":55},"line":{"color":"rgba(31,119,180,1)","width":55},"frame":null}],"highlight":{"on":"plotly_click","persistent":false,"dynamic":false,"selectize":false,"opacityDim":0.2,"selected":{"opacity":1},"debounce":0},"shinyEvents":["plotly_hover","plotly_click","plotly_selected","plotly_relayout","plotly_brushed","plotly_brushing","plotly_clickannotation","plotly_doubleclick","plotly_deselect","plotly_afterplot","plotly_sunburstclick"],"base_url":"https://plot.ly"},"evals":[],"jsHooks":[]}</script>
```

…

Практика и реальный пример

Проведём PCA уже на полных данных pima

# Загрузим библиотеки
#install.packages('devtools')
#remotes::install_github("vqv/ggbiplot")
#install.packages('ggbiplot')

library(devtools)
## Загрузка требуемого пакета: usethis
library(FactoMineR)
library(ggbiplot)
## Загрузка требуемого пакета: plyr
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
## 
## Присоединяю пакет: 'plyr'
## Следующие объекты скрыты от 'package:plotly':
## 
##     arrange, mutate, rename, summarise
## Следующий объект скрыт от 'package:ggpubr':
## 
##     mutate
## Следующие объекты скрыты от 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## Загрузка требуемого пакета: scales
## Загрузка требуемого пакета: grid
## 
## Присоединяю пакет: 'ggbiplot'
## Следующий объект скрыт от 'package:ggfortify':
## 
##     ggbiplot

Делаем PCA:

pima_full.pca <- prcomp(pima_clear, 
                        scale = T) # Не забываем про стандартизацию!

Оценим результат.

summary(pima_full.pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6     PC7
## Standard deviation     1.5999 1.2477 1.0949 0.9776 0.84863 0.63358 0.55774
## Proportion of Variance 0.3199 0.1946 0.1499 0.1195 0.09002 0.05018 0.03888
## Cumulative Proportion  0.3199 0.5145 0.6644 0.7839 0.87387 0.92404 0.96293
##                            PC8
## Standard deviation     0.54458
## Proportion of Variance 0.03707
## Cumulative Proportion  1.00000

Смотрим на “Cumulative Proportion”. У нас в данных первые 4 главные компоненты объясняют 74% вариации данных. Посмотрим это на графике:

factoextra::fviz_eig(pima_full.pca, 
         addlabels = T, 
         ylim = c(0, 40))

На самом деле, это не слишком хороший результат, т.к. следующая конвенционально важная отметка в 90% достигается уже только на PC7. Бывают данные, которые не слишком хороши для PCA, но нельзя сказать, что у нас всё ужасно. Первые две компоненты объясняют 50% дисперсии. Вокруг них и будет сосредоточен основной анализ.

Анализ переменных по PCA

Мы можем начать анализировать, как наши переменные связаны с PC1 и PC2. Посмотрим на график ниже:

factoextra::fviz_pca_var(pima_full.pca, col.var = "contrib")

Подсказка для интерпретации графика здесь.

  • Стрелки - средние значения переменных для PC1 (Dim1) и PC2 (Dim2). В скобках указаны проценты объяснённой дисперсии каждой из двух компонент. На каждую последующую PC всегда приходится всё меньше и меньше разброса в данных.
  • Цвет и близость к кругу - насколько та или иная переменная вносит вклад в анализируемые главные компоненты
  • Направление - относительная мера близости переменных. Если стрелки расходятся в прямо-противоположные стороны, то переменные отрицательно скоррелированы внутри представленных главных компонент.

В данных мы видим три группы переменных: * age, pregnant. * mass, triceps, pedigree * остальные

Теперь оценим PC1 и PC2 в целом. Что "схватили" первая и вторая компоненты? какой разброс мы видим по осям x и y?

…

Мы также можем отдельно посмотреть на, например, топ 3 самых важных переменных с т.зр. их вариации в PC1 и PC2:

factoextra::fviz_pca_var(pima_full.pca, 
             select.var = list(contrib = 3), # Задаём число компонент здесь 
             col.var = "contrib")

У самих по себе главных компонент есть одна очень большая проблема с точки зрения их анализа. Какая?

Посмотрим из чего состоят 1, 2 и 3 главные компоненты:

factoextra::fviz_contrib(pima_full.pca, choice = "var", axes = 1, top = 24) # 1

factoextra::fviz_contrib(pima_full.pca, choice = "var", axes = 2, top = 24) # 2

factoextra::fviz_contrib(pima_full.pca, choice = "var", axes = 3, top = 24) # 3

Анализ наблюдений по PCA

Помимо переменных мы можем анализировать также и наблюдения, искать в них кластеры и корреляцию с переменными в целом. Для этого используется biplot. Bi потому, что на нем одновременно изображены и точки, и переменные

# Загрузим библиотеку
library(ggbiplot)

Сделаем biplot:

ggbiplot(pima_full.pca, 
         scale=0, alpha = 0.1) + 
  theme_minimal()

Более осмысленным biplot становится при использовании кластерных методов, с помощью которых мы можем разделить наблюдения на группы. Посмотрим, наблюдается ли разница между группами по diabetes:

# Сделаем корректные данные для группировки по diabetes.
pima_clear_with_ch <- pima %>% 
  filter(glucose != 0 & pressure != 0 & triceps != 0 & insulin != 0 & mass != 0 & age != 0 )

# Визуализируем с группировкой по diabetes (для этого переменную нужно сделать фактором)
ggbiplot(pima_full.pca, 
         scale=0, 
         groups = as.factor(pima_clear_with_ch$diabetes), 
         ellipse = T,
         alpha = 0.2) +
  theme_minimal()

А что с возрастными группами:

ggbiplot(pima_full.pca, 
         scale=0, 
         groups = as.factor(pima_clear_with_ch$age_group), 
         ellipse = T,
         alpha = 0.2) +
  theme_minimal()

Для дальнейшего ознакомления с PCA я рекомендую посмотреть следующие туториалы: * https://bioinfo4all.wordpress.com/2021/01/31/tutorial-6-how-to-do-principal-component-analysis-pca-in-r/ * https://towardsdatascience.com/principal-component-analysis-pca-101-using-r-361f4c53a9ff * https://juliasilge.com/blog/best-hip-hop/ * https://juliasilge.com/blog/cocktail-recipes-umap/ * http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/118-principal-component-analysis-in-r-prcomp-vs-princomp/

…

Завершение. А если графики используются не для эксплораторного анализа, а для презентации результатов?

В визуализации данных есть своя теория и свои исследования воприятия разных графиков (например что-то вы можете почитать на data-to-viz в блоге или у Анастасии Кузнецовой).

Когда график для вас - способ представить результат вашим коллегам или широкой публике, его стоит делать исходя из несколько других соображений, чем когда вы делаете это при эксплораторном анализе. По этой причине, при подготовке графика следует учесть следующие принципы:

  1. Оцените насколько хорошо ваши данные подходят типу графика.

  2. Делайте фокус на чем-то одном: общих паттернах или деталях. Исходя из этого стоит выбирать тип графика.

  3. Агрегируйте большие объемы данных при визуализации.

  4. Правильно выбирайте палетки. Учитывайте то, как их могут читать дальтоники или, например, акцентируют ли цвета внимание читателя на том, что вам нужно?

  5. Не делайте график впечатлительным без необходимости. Эффектность в простоте.

  6. Убедитесь, что график соответствует интуитивным конвенциям восприятия (у него не перевёрнуты оси, они не обрезаны (но иногда это позволительно), у данных указан источник, все подписи унифицированны и проч.)

  7. При интерпретации результатов помните - график показывает только то, что он показывает. Ни один график не показывает вам причинных эффектов. Только связи. Correlation != Causation.

  8. Помните, что цель любого хорошего графика – рассказать историю (и убедить вас в ней).

…

Я попытался рассказать историю на наших данных на графике ниже. Попробуйте и вы!

plot <- pima %>% 
  mutate(
    age_group = factor(age_group, levels = c("21-30", "31-40", "41-50", "51-60", "60+" )),
    diabetes = case_when(
      diabetes == 'pos' ~ "Diabet-Positive",
      diabetes == 'neg' ~ "Diabet-Negative"
    ) # Переводим в фактор, для корректной последовательности категорий
    ) %>% 
  filter(insulin != 0 & glucose != 0) %>% 
  ggplot(aes(x=insulin, y=glucose, color = age_group)) + 
  geom_point(size = 3, alpha = 0.8) + 
  facet_grid(. ~ diabetes) +
  scale_color_brewer(palette = 'OrRd') +
  ggtitle('Indian Women Diabets') +
  labs(y = 'Plasma glucose concentration (log10)', x = '2-Hours serum insulin (mu U\\ml) (log10)') + 
  guides(color = guide_legend(title = 'Age Groups')) +
  scale_x_log10() + scale_y_log10() +
  theme_minimal() 

plot

Полезные ссылки по теме визуализации данных